home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / VideoToolbox 96.06.15 / VideoToolboxSources / NoisePdfFill.c < prev    next >
Text File  |  1996-04-03  |  17KB  |  420 lines
  1. /*
  2. NoisePdfFill.c
  3.  
  4. error=NoisePdfFill(world,&rect,pdfKind,&mean,&sd,min,max);
  5. error=NoisePdfAdd(world,&rect,pdfKind,&mean,&sd,min,max);
  6.  
  7. Generate good noise images quickly. You supply a GWorld (or a window) and a rect
  8. indicating how much of it to fill with noise. "world" may be a pointer to a
  9. GWorld, a color CWindow, or, with casting, a plain-old Window. NoisePdfFill fills
  10. the pixels within rect with noise; NoisePdfAdd adds noise to them. Each pixel
  11. will get an independent integer sample from a specified probability density
  12. function (pdf): kBinaryPdf, kUniformPdf, or kGaussianPdf. (A fourth pdf,
  13. kBinomialPdf, is not recommended.) You specify the pdf kind, mean, and
  14. "unclipped" standard deviation. The last two arguments, "min" and "max", are used
  15. to clip, i.e. restrict the range of the pdf (e.g. to preclude overflow in your
  16. image pixels). In the cases of kUniformPdf, kBinomialPdf, and kGaussianPdf,
  17. samples falling outside the range [min,max] are discarded and replaced by
  18. independent new samples. In the case of kBinary, each sample has one of only two
  19. values, with equal probability; these two values are initially set to mean±sd, to
  20. yield the desired mean and sd, but they are adjusted, if necessary, to bring them
  21. within the range [min,max]. Upon return, the mean and sd arguments (which are
  22. passed by reference, i.e. by address, not by value) return the mean and standard
  23. deviation of the random generator of the samples (i.e. the clipped pdf).
  24.  
  25. The requested standard deviation, sd, must be greater than or equal to zero. The
  26. special case of zero standard deviation is recognized and performed quickly,
  27. producing a uniform image with specified mean. A request for an unachievably low,
  28. but nonzero, sd will be served by providing the lowest achievable sd (given the
  29. specified min and max), on the premis that the user cares most about log sd, so
  30. that the least possible nonzero sd will be a better approximation to an
  31. unachievably low sd than would zero sd. Supplying an sd of INF will produce
  32. samples with values of min and max if the pdf kind is kBinaryPdf, and will
  33. produce a uniform distribution over the range [min,max] if the pdf is
  34. kUniformPdf, kGaussianPdf, or kBinomialPdf.
  35.  
  36. "max" must be greater than or equal to "mean", which must be greater than or
  37. equal to "min". Normally, image pixels are considered positive and you'll want to
  38. make "min" greater than or equal to zero, but this is not enforced by
  39. NoisePdfFill; e.g. you can add zero-mean noise to a positive image. 
  40.  
  41. Calling NoisePdfAdd with a "min"<0 should work perfectly (though it hasn't been
  42. tested), but may be significantly slower than when "min"≥0. NoisePdfAdd() uses
  43. NoisePdfFill() to make the noise in a temporary new GWorld and then calls
  44. CopyWindows to add the noise to your image. If "min"≥0 then NoisePdfAdd() uses a
  45. fast parallel addition that assumes there will be no overflow in each pixel's
  46. arithmetic, so the user should take care to specify a "max" that will indeed
  47. prevent overflow, since any overflow would affect a neighboring pixel. If "min"<0
  48. then the unsigned addition of "negative" with positive pixels to produce positive
  49. pixels will "overflow" in the normal course of computing the sum using unsigned
  50. arithmetic. The bits of the pixel itself will be correct, but the overflow would
  51. affect a neighboring pixel if parallel addition were used, so the additions are
  52. done independently, one pixel at a time, which takes longer.
  53.  
  54. kBinaryPdf, kUniformPdf, and kGaussianPdf should all be useful. We use
  55. kGaussianPdf in most of our noise-masking work, but for dynamic noise we
  56. sometimes use kBinaryPdf in order to achieve more noise power, assuming that the
  57. observer will visually integrate several of the dynamic frames, blurring the
  58. distinction between pdf's. Presumably the three kinds of noise would be visually
  59. indistinguishable if shown with independent 15 ms frames and equal sd.
  60.  
  61. kBinomialPdf works fine, but is not recommended. The kBinomialPdf code uses the
  62. fact that the number of heads in 24 coin tosses is approximately gaussian. This
  63. code was originally written as a quick way to make approximately gaussian noise.
  64. However, kBinomialPdf is probably much slower than the subsequently written
  65. kGaussianPdf, which provides samples that are accurately gaussian.
  66.  
  67. void GetGoodNoisePdfBounds(pdfKind,mean,clippedSd,&unclippedSd,&min,&max);
  68.  
  69. GetGoodNoisePdfBounds is a utility routine that helps you choose good parameters
  70. with which to call NoisePdfFill or NoisePdfAdd. Be aware that, whereas
  71. NoisePdfFill and NoisePdfAdd are general-purpose, the "good" in the name
  72. "GetGoodNoisePdfBounds" indicates that its recommendations are not necessarily
  73. best. It incorporates several reasonable but somewhat arbitrary assumptions, e.g.
  74. that you'd like to clip your gaussian at ±2 sd, and helps you to choose a good
  75. set of parameters to obtain that result. GetGoodNoisePdfBounds will not be
  76. appropriate to all uses of NoisePdfFill and NoisePdfAdd; you may want to write
  77. your own. Note: GetGoodNoisePdfBounds is still under development; I haven't yet
  78. added its prototype to VideoToolbox.h. Beware that if the ±2*sd bounds are 
  79. clipped by the supplied min,max bounds then the returned unclippedSd will 
  80. be wrong. I wanted to add code to compute the exact answer, when it exists, 
  81. but the integrals for the gaussian case were sufficiently difficult that I 
  82. put it aside for some future date.
  83.  
  84. The hardest part of writing the code in this file was figuring out how to
  85. separate the general-purpose code (i.e. NoisePdfFill) from the pragmatic but
  86. somewhat arbitrary assumptions now encapsulated in GetGoodNoisePdfBounds.
  87.  
  88. HISTORY:
  89. 3/20/95 dgp wrote it, based on my MakeTextInNoiseWorld.c
  90. 4/8/95 dgp removed any assumption that world is a GWorldPtr; it can also
  91. be a CWindowPtr or a WindowPtr.
  92. 4/18/95 dgp fixed error in NoisePdfAdd: mean & sd were swapped.
  93. */
  94.  
  95. #include "VideoToolbox.h"
  96. OSErr NoisePdfFill(GWorldPtr world,Rect *rectPtr
  97.     ,int pdfKind,double *meanPtr,double *sdPtr,int min,int max);
  98. OSErr NoisePdfAdd(GWorldPtr world,Rect *rectPtr
  99.     ,int pdfKind,double *meanPtr,double *sdPtr,int min,int max);
  100. void GetGoodNoisePdfBounds(int pdfKind,double mean,double clippedSd
  101.     ,double *unclippedSdPtr,int *minPtr,int *maxPtr);
  102.  
  103. //enum{kBinaryPdf=0,kBinomialPdf,kGaussianPdf,kUniformPdf}; // in VideoToolbox.h
  104.  
  105. #undef MAX
  106. #undef MIN
  107. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  108. #define MIN(a,b) ((a) < (b) ? (a) : (b))
  109. static Boolean diagnostics=0;
  110.  
  111. OSErr NoisePdfFill(GWorldPtr world,Rect *rectPtr
  112.     ,int pdfKind,double *meanPtr,double *sdPtr,int min,int max)
  113. {
  114.     CGrafPtr oldPort;
  115.     GDHandle oldDevice;
  116.     int error=0;
  117.     OSErr osErr=0;
  118.     long n;
  119.     int i,j,x,y;
  120.     long low,high,sample;
  121.     double actualMean,actualSd;    // mean and sd of process that generated the image
  122.     long actualMin,actualMax;
  123.     int binomialSamples;
  124.     double binomialFactor;
  125.     long extraSamples;
  126.     short *distribution;
  127.     static short **distributionHandle=NULL;
  128.     static long distributionSize;
  129.     static double lowSave,highSave,meanSave,unclippedSdSave,actualSdSave;
  130.     double meanSquare,a,aa;
  131.     RGBColor rgb;
  132.     GDHandle device;
  133.     int pixelSize;
  134.     static const RGBColor blackRGB={0,0,0},whiteRGB={0xFFFF,0xFFFF,0xFFFF};
  135.     unsigned long pix[2048];
  136.     register long *pL;
  137.     
  138.     // Make noise. Takes 300 ms
  139.     assert(StackSpace()>4000);
  140.     assert(world!=NULL && sdPtr!=NULL && meanPtr!=NULL);
  141.     assert(!IsNan(*meanPtr) && !IsNan(*sdPtr));
  142.     assert(*sdPtr>=0.0);
  143.     assert(min<=*meanPtr && *meanPtr<=max);
  144.     high=max;
  145.     low=min;
  146.     if(IsInf(*sdPtr) && pdfKind!=kBinaryPdf)pdfKind=kUniformPdf;
  147.     switch(pdfKind){
  148.     case kBinaryPdf:
  149.         // high and low are values corresponding to heads and tails.
  150.         // The high-low difference is more important than their individual values.
  151.         if(!IsInf(*sdPtr)){
  152.             a=2*(*sdPtr);
  153.             low=MAX(round(*meanPtr-a/2),min);
  154.             high=MIN(round(low+a),max);
  155.             low=MAX(round(high-a),min);
  156.             if(low==high && *sdPtr>0){
  157.                 if(high<max)high++;
  158.                 else if(low>min)low--;
  159.             }
  160.         }
  161.         break;
  162.     case kUniformPdf:
  163.         if(!IsInf(*sdPtr)){
  164.             a=2*sqrt(3)*(*sdPtr)-1;    // integer range of pdf for specified sd
  165.             low=MAX(round(*meanPtr-a/2),min);
  166.             high=MIN(round(low+a),max);
  167.             low=MAX(round(high-a),min);
  168.             if(low==high && *sdPtr>0){
  169.                 if(high<max)high++;
  170.                 else if(low>min)low--;
  171.             }
  172.         }
  173.         break;
  174.     default:
  175.         break;
  176.     }
  177.     GetGWorld(&oldPort,&oldDevice);
  178.     device=GetWindowDevice((WindowPtr)world);
  179.     pixelSize=(**(**device).gdPMap).pixelSize;
  180.     if(low==high){
  181.         SetGWorld(world,device);
  182.         Index2Color(low,&rgb);
  183.         RGBBackColor(&rgb);
  184.         EraseRect(rectPtr);
  185.         SetGWorld(oldPort,oldDevice);
  186.         actualMean=low;
  187.         actualSd=0;
  188.         goto done;
  189.     }
  190.     switch(pdfKind){
  191.     case kBinaryPdf:
  192.         SetGWorld(world,device);
  193.         Index2Color(low,&rgb);
  194.         RGBBackColor(&rgb);
  195.         Index2Color(high,&rgb);
  196.         RGBForeColor(&rgb);
  197.         error=NoiseFill(world,rectPtr,1,1,0); 
  198.         RGBBackColor(&whiteRGB);
  199.         RGBForeColor(&blackRGB);
  200.         SetGWorld(oldPort,oldDevice);
  201.         actualSd=0.5*(high-low);
  202.         actualMean=0.5*(high+low);
  203.         break;
  204.     case kUniformPdf:
  205.         n=1+high-low;
  206.         assert(rectPtr->right-rectPtr->left<=sizeof(pix)/sizeof(*pix));
  207.         for(y=rectPtr->top;y<rectPtr->bottom;y++){
  208.             pL=(long *)pix;
  209.             for(x=rectPtr->right-rectPtr->left-1;x>=0;x--)
  210.                 *pL++=low+nrand(n);
  211.             SetWindowPixelsQuickly((WindowPtr)world,rectPtr->left,y
  212.                 ,pix,rectPtr->right-rectPtr->left);
  213.         }
  214.         assert(n<=sizeof(pix)/sizeof(*pL));
  215.         pL=(long *)pix;
  216.         for(x=low;x<=high;x++)*pL++=x;
  217.         actualMean=MeanL((long *)pix,1+high-low,&actualSd);
  218.         if(diagnostics)printf("unclippedSd %f, actualSd %f, actualSd/unclippedSd %f\n"
  219.             ,*sdPtr,actualSd,actualSd/(*sdPtr));
  220.         break;
  221.     case kGaussianPdf:
  222.         /*
  223.         First we make an ordered list of about twenty thousand integers, with 
  224.         frequency proportional to the desired probability. This table only needs
  225.         to be made once and can be reused indefinitely, until the parameters  
  226.         (mean, sd, min, max) of the distribution change. Each noise pixel is
  227.         a random sample from the table.
  228.         */
  229.         if(distributionHandle==NULL || low!=lowSave || high!=highSave 
  230.             || IsInf(*sdPtr)!=IsInf(unclippedSdSave)
  231.             || (!IsInf(*sdPtr) && fabs(*sdPtr/unclippedSdSave-1)>.01)
  232.             || *meanPtr!=meanSave){
  233.             distributionSize=256*100L;
  234.             if(distributionHandle==NULL){
  235.                 distributionHandle=(short **)TempNewHandle(distributionSize*sizeof(*distribution),&osErr);
  236.                 error=osErr;
  237.                 if(distributionHandle==NULL){
  238.                     distributionHandle=(short **)NewHandle(distributionSize*sizeof(*distribution));
  239.                     error=MemError();
  240.                 }else error=0;
  241.                 if(distributionHandle==NULL)PrintfExit("Sorry, couldn't allocate %ld bytes "
  242.                     "for a gaussian table. File “%s” line %d.\n"
  243.                     ,distributionSize*sizeof(*distribution),__FILE__,__LINE__);
  244.             }
  245.             if(fabs((high-*meanPtr)-(highSave-meanSave))>0.5 || high-low!=highSave-lowSave 
  246.                 || IsInf(*sdPtr)!=IsInf(unclippedSdSave)
  247.                 || (!IsInf(*sdPtr) && fabs(*sdPtr/unclippedSdSave-1)>.01)){
  248.                 // takes 300 ms
  249.                 BoundedNormalIntegers(*distributionHandle,distributionSize
  250.                     ,*meanPtr,*sdPtr,low,high);
  251.             }else{
  252.                 j=round(*meanPtr-meanSave);
  253.                 distribution=*distributionHandle;
  254.                 if(j!=0)for(i=0;i<distributionSize;i++)distribution[i]+=j;
  255.                 *meanPtr=meanSave+j;
  256.             }
  257.             lowSave=low;
  258.             highSave=high;
  259.             unclippedSdSave=*sdPtr;
  260.             meanSave=*meanPtr;
  261.             actualMean=MeanW(*distributionHandle,distributionSize,&actualSdSave);
  262.         }
  263.         actualSd=actualSdSave;
  264.         assert(rectPtr->right-rectPtr->left<=sizeof(pix)/sizeof(*pix));
  265.         //HLock((Handle)distributionHandle);
  266.         distribution=*distributionHandle;
  267.         for(y=rectPtr->top;y<rectPtr->bottom;y++){
  268.             pL=(long *)pix;
  269.             for(x=rectPtr->right-rectPtr->left-1;x>=0;x--)
  270.                 *pL++=distribution[nrand(distributionSize)];
  271.             SetWindowPixelsQuickly((WindowPtr)world,rectPtr->left,y
  272.                 ,pix,rectPtr->right-rectPtr->left);
  273.         }
  274.         if(diagnostics)printf("unclippedSd %f, actualSd %f, actualSd/unclippedSd %f\n"
  275.             ,*sdPtr,actualSd,actualSd/(*sdPtr));
  276.         break;
  277.     case kBinomialPdf:
  278.         // Use the number of heads in 24 coin flips as a random variable.
  279.         // It can be computed quickly and is approximately normal.
  280.         // Shift and scale to obtain the specified mean and variance.
  281.         // Discard samples outside the range [min,max].
  282.         // Return mean and sd of the clipped distribution.
  283.         binomialSamples=24;    // should be even (for integer divide by 2)
  284.         binomialFactor=*sdPtr/sqrt(binomialSamples*0.25);
  285.         /*
  286.         The empirical approach of actually measuring the mean and variance
  287.         of our noise generator demands that we produce a
  288.         reasonable number of samples, at least 1000. Any extra samples are later
  289.         discarded.
  290.         */
  291.         assert(rectPtr->right-rectPtr->left<=sizeof(pix)/sizeof(*pix));
  292.         extraSamples=ceil(1000.0/(rectPtr->bottom-rectPtr->top))
  293.             -(rectPtr->right-rectPtr->left);
  294.         if(extraSamples<0)extraSamples=0;
  295.         a=aa=0.0;
  296.         n=0;
  297.         for(y=rectPtr->top; y<rectPtr->bottom; y++){
  298.             pL=(long *)pix;
  299.             for(x=rectPtr->right-rectPtr->left-1; x>=-extraSamples; x--){
  300.                 do{
  301.                     sample=round(*meanPtr+binomialFactor
  302.                         *(BinomialSampleQuickly(binomialSamples)-binomialSamples/2));
  303.                 }while(sample<low || sample>high);
  304.                 a+=sample;
  305.                 aa+=(long)sample*(long)sample;
  306.                 n++;
  307.                 if(x>=0)*pL++=sample;
  308.             }
  309.             SetWindowPixelsQuickly((WindowPtr)world,rectPtr->left,y
  310.                 ,pix,rectPtr->right-rectPtr->left);
  311.         }
  312.         actualMean=a/n;
  313.         actualSd=sqrt(aa/n-actualMean*actualMean);
  314.         if(diagnostics)printf("unclippedSd %f, actualSd %f, actualSd/unclippedSd %f, n %ld\n"
  315.             ,*sdPtr,actualSd,actualSd/(*sdPtr),n);
  316.         break;
  317.     default:
  318.         PrintfExit("%s,%d: No such pdf %d\n",__FILE__,__LINE__,(int)pdfKind);
  319.     }
  320.     if(diagnostics){
  321.         ImageStatistics(world,rectPtr,&actualMin,&actualMax,&actualMean,&meanSquare);
  322.         if(actualMax>max || actualMin<min)printf("\nIMAGE OVERFLOW ERROR:\n");
  323.         printf("noise actualMin %ld≥%d, actualMax %ld≤%d, actualMean %.1f≈%.1f, rms %.1f≈%.1f≈%.1f\n"
  324.             ,actualMin,(int)min
  325.             ,actualMax,(int)max
  326.             ,actualMean,*meanPtr
  327.             ,sqrt(meanSquare),actualSd,*sdPtr);
  328.         assert(actualMax<=max && actualMin>=min);
  329.     }
  330. done:
  331.     *sdPtr=actualSd;
  332.     *meanPtr=actualMean;
  333.     return error;
  334. }
  335.  
  336. OSErr NoisePdfAdd(GWorldPtr world,Rect *rectPtr
  337.     ,int pdfKind,double *meanPtr,double *sdPtr,int min,int max)
  338. {
  339.     GWorldPtr noiseWorld;
  340.     GDHandle device;
  341.     int error=0,pixelSize;
  342.     double meanSquare,actualMean;
  343.     long actualMin,actualMax;
  344.     ColorTable **ct;
  345.     long copyMode;
  346.  
  347.     assert(world!=NULL);
  348.     device=GetWindowDevice((WindowPtr)world);
  349.     pixelSize=(**(**device).gdPMap).pixelSize;
  350.     ct=GetCTable(pixelSize+32);    // gray ramp
  351.     assert(ct!=NULL);
  352.     error=NewGWorld(&noiseWorld,pixelSize,rectPtr,ct,0,0);
  353.     if(error)error=NewGWorld(&noiseWorld,pixelSize,rectPtr,ct,0,useTempMem);
  354.     if(error)return error;
  355.     error=NoisePdfFill(noiseWorld,rectPtr,pdfKind,meanPtr,sdPtr,min,max);
  356.     if(error)goto done;
  357.     if(min>=0)copyMode=addOverParallelLiterally;    // no overflow
  358.     else copyMode=addOver;                            // "sign" bit may overflow
  359.     error=CopyWindows(noiseWorld,world,rectPtr,rectPtr,copyMode,NULL);
  360.     if(error)PrintfExit("%s line %d: CopyWindows error %d\n",__FILE__,__LINE__,error);
  361.     if(diagnostics){
  362.         ImageStatistics(world,&world->portRect,&actualMin,&actualMax,&actualMean,&meanSquare);
  363.         if(actualMax>max || actualMin<min || diagnostics){
  364.             printf("signal+noise actualMin %ld, actualMax %ld, actualMean %.1f, actualSd %.1f\n"
  365.                 ,actualMin,actualMax,actualMean,sqrt(meanSquare-actualMean*actualMean));
  366.         }
  367.     }
  368. done:
  369.     DisposeGWorld(noiseWorld);
  370.     return error;
  371. }
  372.  
  373. void GetGoodNoisePdfBounds(int pdfKind,double mean,double clippedSd
  374.     ,double *unclippedSdPtr,int *minPtr,int *maxPtr)
  375. {
  376.     double dMax,dMin;    // nominal range of difference from mean
  377.     double clippedSd1;    // value of clippedSd that would be produced if unclippedSd were 1
  378.     int low,high;        // range of difference from mean
  379.     
  380.     assert(unclippedSdPtr!=NULL && minPtr!=NULL && maxPtr!=NULL);
  381.     assert(clippedSd>=0);
  382.     assert(*maxPtr>=*minPtr);
  383.     
  384.     switch(pdfKind){
  385.     case kBinaryPdf:
  386.         dMax=clippedSd;
  387.         dMin=-dMax;
  388.         clippedSd1=1;
  389.         break;
  390.     case kUniformPdf:
  391.         dMax=sqrt(3)*clippedSd;
  392.         if(dMax>=1)dMax-=0.5;
  393.         dMin=-dMax;
  394.         clippedSd1=1;
  395.         break;
  396.     case kBinomialPdf:
  397.         dMax=2*clippedSd;
  398.         dMin=-dMax;
  399.         clippedSd1=0.85;
  400.         break;
  401.     case kGaussianPdf:
  402.         dMax=2*clippedSd;
  403.         dMin=-dMax;
  404.         clippedSd1=0.737;
  405.         break;
  406.     }
  407.     // high and low are the extremes of the pdf.
  408.     // The high-low difference is more important than the absolute value of either.
  409.     assert(dMax>=dMin);
  410.     low=MAX(round(mean+dMin),*minPtr);
  411.     high=MIN(round(low+(dMax-dMin)),*maxPtr);
  412.     low=MAX(round(high-(dMax-dMin)),*minPtr);
  413.     if(low==high && clippedSd>0){
  414.         if(high<*maxPtr)high++;
  415.         else if(low>*minPtr)low--;
  416.     }
  417.     *minPtr=low;
  418.     *maxPtr=high;
  419.     *unclippedSdPtr=clippedSd/clippedSd1;
  420. }